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Abstract. Due to Fortuin and Kastelyin the q state Potts model has a representation 
as a sum over random graphs, generalizing the Potts model to arbitrary q is based 
on this representation. A key element of the Random Cluster representation is the 
combinatorial factor r5(C, E), which is the number of ways to form C distinct clusters, 
consisting of totally E edges. We have devised a method to calculate rg;(C, E) from 
Monte Carlo simulations. 
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1. Introduction 

The Potts model P] is one of the most studied models in statistical physics. The 
traditional representation of the model is in terms of the Hamiltonian 



where the spins o"i are integer values o"j G the sum {i,j) is over nearest 

neighbours. The g is a parameter of the model. The model is typically defined on 
a regular lattice in d dimensions, but can in general be defined on any graph. 

For d > 2 the model sustains a order-disorder transition , in d = 2 the critical 
coupling is Pc = ln(l + ^/q). For [3 > (3c the g-fold permutation symmetry of Eq. [T]is 
broken, and one of the q different groundstates has been singled out. For q = 2 the model 
is the familiar Ising model, which has a second order transition, but with increasing q 
the excited states have relatively more entropy and for q > qc the transition is first 
order. For d = 2 the phase transition changes order at gc = 4 El , for c? = 3 the exact 
value is not known, but the most recent estimate based on Monte Carlo simulations is 
qc ^ 2.35|1I. 

The Hamiltonian Eq. ^ is only defined for integer q, however due to an elegant 
transformation by Fortuin and Kastelyn (KF) the partition function of the q state Potts 
model can be written as a correlated percolation problem, the socalled Random Cluster 
(RC) modeljS]. In the RC representation q enters as an ordinary variable, and can 
attain any scalar value. Apart from extrapolation/interpolation from integer q results, 
all (numerical) studies of the noninteger q properties of the Potts model are based 
on the RC representation, this also applies to the current paper. Properties of the 
Potts model with noninteger q have been extensively studied using transfer matrixjH] 
techniques. Recently also MC simulations have been used. The latter come in two 
categories; either a technique is based on the RC measure to simulate directly at an 
arbitrary gj^EllH], or alternatively the results are reweighted to arbitrary q after the 
simulation is complete [3 Ej. 

The rest of this paper is organized as follows: In section |21 we introduce some 
key elements of graph theory, and how concepts from graph theory can be applied 
in statistical physics; in particular to the Potts model. In section El we introduce 
and describe an algorithm which can be used to "reweight" Potts model simulations 
to arbitrary q. Section HI is devoted to results, both to show the correctness of the 
approach and also to study real q properties which are not easily studied by ordinary 
MC simulations. 

2. Graph theory and the Potts model 

An (undirected) graph ^ is a collection of vertices V{Q), along with a set of edges E{Q) 
connecting the vertices jTHI. A subgraph Q' E Q is a. collection of vertices and edges such 




(1) 
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that V{Q') G V{Q) and E{Q') G E{Q). The rank of a graph is denoted by r{Q) and 
given by 

r{G) = \Vm-C{G), (2) 

where |V^(^)| is the number of vertices and C{Q) is the number of connected components. 
Observe that also isolated single vertices constitute connected components when 
evaluating the rank of a graph. Fig.^shows a simple graph and illustrates the necessary 
concepts. From now on we will use the symbols Eg, Cg and \/g to denote the number of 
edges, clusters and vertices in a graph when there is no ambiguity we will omit the 
index Q. 



Figure 1. The sites of 5 x 5 lattice, and links connecting some of the sites. Together 
these sites and links constitute a graph. This particular graph has Vg = 25, Eg = 20, 
six connected components {Cg = 6) and a rank r{Q) = Vg — Cg = 19. 

By assigning scalar properties to sites and bonds one can define different graph 
polynomials. One of the most general graph polynomials is the Tutte or Di-Chromatic 
polynomal Tg {x, y) [TTl [T2j : 

T,(x,y)= {x-if'''-'''Hy-if-''''- (3) 

EeE{g) 

The sum in Eq.Elis over all edge configurations of the graph Q (i.e. spanning subgraphs). 
Here a; is a scalar property assigned to the vertex set, and y a property assigned to the 
edges; as indicated in Eq.Elwe will only consider the situation of spatially constant y, but 
the general definition of the Tutte polynomial allows for a set {y} of edge properties. 
Many other polynomials can be found as suitably rescaled evaluations of the Tutte 
polvnomial|13j: 

Roip) =(l-P)'"''^V"%(l>r^) (4) 
PgiQ) =(-ir(^)g%(l-g,0) (5) 

Zg{q,v) = qv'^-\v + ir^Tg(^,v+l], v = = e^' - 1. (6) 
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Rg{p) is the reliability polynomial, closely related to the (bond) percolation problem. 
Pg{q) is the chromatic polynomial, and denotes the number of ways the vertices in Q can 
be colorized with q different colors, so that no adjacent vertices share the same color. 
The chromatic polynomial coincides with the T — limit of the partition function of 
the antzferromagnetic Potts model. Finally Z{q,v) is the partition function of the q 
state Potts model. Observe the quantity v in Eq. IHl in this context this is the most 
convenient temperature variable. 

The FK transformation is the key to identify Z{q,v) with the Tutte Polynomial|Sj. 
The actual transformtion is in terms of the complete partition function, hence it is not 
possible to identify a spin state with a corresponding RC state uniquely, see however 
Ref. J3] for an exposition in terms of a mixed bond-spin model which elucidates the 
connection. Zrc(p, q) is a function of two variables: a probability p to occupy an edge, 
and a q, where Ing resembles a cluster entropy. The RC partition function is built up as 
follows: (1) each configuration E'(Q) of edges gets a "Boltzmann" -weight —p)^~^\ 
(2) the weight is multiplied by an entropic factor g^', (3) all configurations E'(Q) are 
summed over. This finally gives the RC partition function 

V E 

ZKciq,p) = E - P)'"''^' = E E E)/(l - pf'^' gS (7) 

E'(g)eE{g) c=i e'=o 

V ' 

ac(p) 

The p in Eq. [3 is the probability to occupy an edge, for the RC model this is an 
arbitrary number, however to make contact with the g-state Potts model at coupling j3, 
we must have p = 1 — e~^'^. As indicated in Eq. [3 the partition function can bee seen 
as a polynomial in q, with p dependant coeffiecients. In section 14.41 we will use this to 
determine the zeroes of the partition function in the complex q plane. 

Using the combinatorial factor rg(C, E) to express the sum is the key element in 
Eq. [7| This factor is simply the number of ways to form C connected components with 
E, on the underlying graph Q. This is a purely combinatorial/geometric property which 
can in principle be calculated without any reference to a particular model of statistical 
physics. On the other hand all physical properties are contained in Tg{C, E). Eq. EJalso 
highlights that the Potts model has a common structure independent of q, even though 
the physical properties vary significantly with q. In addition to facilitating the study 
of the Potts model for arbitrary q, the FK representation also serves as the theoretical 
underpinning of the Swendsen-Wang algorithm for spin models [T^ IT3]. 

An important topic in computer science is a formal demarcation of tractable and 
intractable problems. The socalled #P complete problems are counting problems which 
are essentially intractable. Obtaining the partition function of (discrete) system belong 
to this category m E]- Due to this intractability good approximative techniques is 
essential; the Monte Carlo technique is one such approach. Also in computer science 
the use of Monte Carlo techniques to approach NP and ^P complete problems, has 
been popular, see eg. [T7j. Computer scientists Jerrum and Sinclair have devised 
efficient Monte Carlo algorithms (FPRAS) to determine the partition functions of both 
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2D monomer-dimer system, and the 2D Ising modei p^ ll9j. Hence the study of the RC 
and related problems is of interest to scientist from widely different fields. 

3. Algorithm 

The probability -P(e) to find a system in a state with energy e is proportional to g{e)e~^'^, 
where g{e) is the density of states at energy e. That P(e) can be written in this manner 
is the foundation of ordinary e — (3 reweigtingpPl. In the formulation Eq. [7| (p, E) and 
(g, C) are "conjugate" variable pairs; alas Tg{Q., E) can be used to reweight to arbitrary 
q and p; from now on we will mostly use (3 in the text, but it should be understood that 
the relation p = 1 — e~^'^ applies throughout. In the remainder of this section we will 
present an algorithm to estimate rg(C, E) from simulations at different p and q. An 
algorithm based on the same principle was presented by Weigel et. al. in Ref. jH], and 
just recently Hartmann has presented an algorithm based on only (g, C) reweighting|3]. 

The algorithm presented here is general, and will apply to any graph. However for 
ease of notation we have specialized to a two dimensional square lattice with a total 
of N = L X L sites, and 2N edges. The Gibbs probability to find any state with C 
components and E edges is given by|13j: 

r.(C.E)pV--V , 

Zg{q,l3) 

To estimate rg(C, E) we need to generate states distributed according to Eq. |H1 We 
have done this by using the S wendsen- Wang [T3j algorithm on the q state Potts model, 
with integer q. However, one could equally well have used an algorithm generating RC 
states directly[71|H], or alternatively a combination. During the simulation at /i = {q,P) 
a histogram /i^(C, E) is collected. From the histogram /i^q(C,E) we can in principle 
estimate rg(C, E) from Eq. |H1 

f,„(C, E) = E)p-^g-(^^-^)g-S (9) 

where is an (undetermined) normalization constant. rg(C, E) is independent of 
fi, however the estimator in Eq. El has been given index fiQ to indicate that it is based 
on results sampled at these couplings. The estimator Eq. IHlis formally correct, but only 
applicable in a narrow range around the mean values (C)^^ and (E)^;,. By combining 
results obtained at different f3 and q we can get an estimate for Tg{C, E) which is valid 
for a wide range of C and E values. A series of N histograms obtained at couplings 
fii, fi2, ■ ■ ■ , fJ'N can be combined as 

TV 

fg(C,E) = ^^,(C,E)-f,,(C,E), (10) 

i 

where the weight factor Wi{C, E) is given by 

w.{C,E) = J^'fl^^ (11) 
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The normalization constants ^j, i > 1 are determined by maximizing, the (weighted) 
overlap between (the logarithm of) the estimates F^. (C,E). Mathematically this 
amounts to minimizing 

^' = EEE^M,:(C,E)V(C,E)x 

i j>i C,E 

( (e. + In V,(C, E) - Clng,) - (e, + \nh^^{C, E) - C Ing,) (12) 

with ^1 initially fixed at an arbitrary value. The final normalization constant is 
determined by the overall normalization 

5^rg(C,E') = 2^ (13) 

C,E' 

The actual solution of the minimization problem Eq. is found as the solution of a 
system of linear equations. As long as all the histograms h^^{C, E) have finite overlap 
with at least one other histogram h^.{C, E) the solution will be found. The method is 
a generalization of an existing algorithm to determine the density of states 5'(e) [ZH I221- 
Due to the nonlinear nature of the algorithm it is difficult to calculate errors by 
the use of error-propagation. Furthermore the estimation of rg(C, E) is quite time 
consuming, hence computer-intensive methods like Jack-Knife and Bootstrap are not 
very suitable. In the current paper error estimates have been calculated by comparing 
the results from independent simulations. 



4. Results 



Basic thermodynamic results 

In this section we will show how simulations performed at one value qi can be reweighted 
to another q2 7^ qi- Fig. El shows thermodynamics for a g = 4 Potts model. The solid 
line is data obtained at g = 4, and the symbols represent results reweighted from g = 2 
and g = 8 respectively. 



4-.2. The average trajectory in clusters - links space 



In the Random Cluster formalism the state of the system is given by C and E, and it is 
interesting to see how these quantities evolve when the Potts model parameters (3 and 
g are varied. For a fixed value of E the conditional probability P(C|E) is independent 
of /3; hence we can easily plot the mean path the system will follow in (C, E) space. In 
Fig. El we show the conditional mean 

ScC-r(c.E),c 

along with the contours of P(C, E) at the critical coupling, for two different values of 
g. As we can see from Fig. IHl the g behaviour of C and E can conveniently be divided 
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Figure 2. This figure shows from top to bottom free energy, internal energy and 
spesific heat for the q — 4: Potts model, system size is 16 x 16. The solid line shows 
result obtained from a simulation at g = 4, the symbols show results "reweighted" 
from q = 2 and q — 8 respectively. 



in three regions: (1) a low T region where (C|E) ^ 1 quite independent of E, a high T 
region where (C|E) > — E and an intermediate region containg the critical point. It 
is only in the intermediate region there is significant q dependence. 

The contours in Fig. El show the probability density P(C, E) at the critical point, 
for q = 2 and q = 8. The q "reweighting" has similar limitations as ordinary thermal 
reweighting, the statistics is best at the original q value, and can not be extended to 
regions of (C, E) space which have not been sampled. As we can from Fig. Elthe overlap 
between the q = 2 and g = 8 results is very small; hence reweighting between these two 
q values would give unreliable results. 

From Fig. iniwe see that the fluctuations are quite assymetric; they are much larger 
along the direction given by the mean path Eq.Elthan orthogonal to it. The conditional 
distribution function 

, I . r(C,E)g^ 

C E = 15 
Ecr(C,E)gC 

is well described by a Gaussian with width asiq)- The width scales with the number of 
sites as iV^/^, hence the relative fluctuations in the number of clusters scales as iV~^/^ 
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Figure 3. Contour plot of the density P(C, E) at the critical point, for q = 2 and 
q = 8 for a 16 X 16 lattice. The dashed lines show (C|E), which corresponds to the 
path followed in C, E space when temperature is varied. 

and consequently the system will follow an increasingly well defined line in (C, E) space 
when the system size increases. Fig. |3] shows the distribution of the cluster density 
c = C/N for a given link density e = E/N, and finite size scaling of the width of this 
distribution, o"e(g) = a^{q)/N. 




Figure 4. The left figure shows the conditinal distribution P{c\e) at e = (1 — 1/(1 + 
y/qfj/L^, i.e. the critical link density, for the 9 = 3 model. The right figure shows the 
width of the distribution P(c|e) as a function of L, all the curves show a decay. 
The curves for q > 3 have been shifted for clarity. 

In the RC model each cluster can be in q different configurations, hence we get an 
additive entropy contribution of In q from every cluster. Consequently we see that for a 
fixed number of links the average number of clusters will increase with q. On the other 
hand larger amount of entropy per cluster, means that for high q entropy will dominate 
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the competetion between internal energy and entropy at a lower number of clusters, 
and consequently at the critical point (C) decreases with increasing q. These points are 
illustrated in Fig. 
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Figure 5. The mean number density of clusters as a function of q, for a fixed density 
of links and at the {q dependant) critical link density. The results in the figure are 
from a 16 X 16 lattice. 



4-3. Evaluation of the Tutte polonymial 

The Tutte polynomial can be defined in terms of a recursive definition which 
immediately leads to a simple and exact algorithm for computation of Tg(x, y). However 
this algorithm has exponential complexity, and is clearly not feasible for anything but 
very small graphs. Due to it's importance in many different areas of mathematics and 
computer science, this has lead to a large effort to find efficient approximate algorithms 
for evaluation of the Tutte polynomial [23j. 

Using the algorithm presented here we can also estimate Tutte polynomials, in 
Fig. ini we show the reliability polynomial and the Chromatic polynomial. With the 
current approach the running time to determine the Tutte polynomial is governed by 
the running time of the MC algorithm, and at least for g < 4 the Swendsen-Wang 
algorithm is rapidly mixing |243. 

When the arguments x,y of the Tutte polynomial move a long way away from the 
values used when samphng, the results become unreliable; consult Eq. IHl to see how 
X and y are related to the parameters q and (3 of the Potts model. In particular for 
a; < 1 and/or y < I the evaluation of T{x,y) is difficult, because in these regions the 
polynomial terms are oscillating and inaccurate coefficients lead to large relative errors. 

4.4- Zeros in the complex q plane 

The formulation of the partition function as a polynomial in q allows for quite easy 
evaluation of the zeros of the partition function in the complex q plane. Properties 
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Figure 6. The reliability polynomial Eq.01and chromatic polynomial Eq.^lfor a 3 x 3 
lattice. The solid lines are exact results from the computer algebra system Maple, and 
the points come from our simulations. The very small system size considered is to 
limit the run-time of Maple 



of the complex q zeroes have been investigated both analytically, and numerically 
According to the Yang-Lee view of critical phenomena the critical point is characterized 
by zeros in the complex f3 plane pinching the real axis. The phase transition in the 
Random Cluster model can be driven by both j3 and q, we should therefor see the same 
pinching of the real q axis. 

The critical coupling is given by jScJ = ln(l + ^/q), alternatively we find that for a 
fixed P the critical q is given by 

= {e^' - if = v\ (16) 

For the current discussion the temperature variable w, first introduced in Eq. IHlwill be 
the most convenient. Plotting the zeros of Z{v, q) we expect the zeros to pinch the 
real q axis close to the qc given by Eq. Fig. [7| shows the distribution of zeros in the 
complex q plane for two different couplings. 

If we denote the zero closest to qc with qdL), we find that qdL) converges towards 
qc with increasing system size. To determine which zero is indeed the "critical" one we 
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Figure 7. The roots in the complex q plane of the partition function Z{v, q) at 
couphngs V — -\/3 (top) and v = (bottom). We observe that the zeroes close in on 
the critical q values of 3 and 5. 

have measured distance d{qi, qc) using both the ordinary metric d2{x,y) = \x — y\ and 
also di{x,y) = |Im(x) — lm{y)\. For v"^ < 3.0 the two methods select the same zero, 
whereas for f ^ > 3.0 different zeros are selected, and the real part of the zero selected 
by d2 jumps about randomly. Fig. |H1 shows finite size scaling plots of the |Im(g)| (as 
determined by using di) for the zero closest to the real q axis. This should scale as 

\lm{q)\ L-^ . (17) 

For q = 2 and g = 3 this gives u ~ 0.992(7) and u ^ 0.863(7) which agree reasonably 
well with the exact values of 1 and 5/6 ~ 0.8333 . . .. For g = 4 we get u ^ 0.77(3), this 
is well above the exact value of 2/3 + logarithmic corrections. If we assume an effective 
exponent for the first order transition at g = 5 we would expect u = 1/2, whereas the 
estimated value is z/ = 0.77(6). 

The reason that the quality of the u estimates detoriate with increasing q is 
probably that the slope of the curve l3c{q) is reduced with increasing q. When the 
transition is driven by q the critical point is approached more and more tangentially. 
It seems reasonable that this makes a precise determination of the critical properties 
progressively more difficult. Furthermore the model has limiting behaviour at g = 4, 
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Figure 8. The plots show |Im((7)| for the zero closest to the real axis, as a function of 
system size L. The value of coincides with Qc. The error bars are generally smaller 
than the symbol size. The solid lines are least squares fits with slope, from top to 
bottom, -1.3(1), -1.29(5), -1.159(9), -1.008(7). 



The zeroes are found using the MPSolve||27j package. To determine the roots of 
Z{v, q) in the complex q plane is an ill-posed problem. Firstly the coefficeints ac{p) (see 
Eq. [Tj) vary over a wide range, secondly finite sampling statistics adds to the problem. 
In particular the states with C — > are typically not sampled at all. For independent 
simulations the pattern of zeroes differs significantly from case to case, however the 
location of the zero q^L) shows much less fluctuations. The results in Fig. |H1 are the 
total of ten independent simulations, and as we see the error bars are very small. 

In a large paper by Alan Sokalj2S] it is shown that the complex q zeros of the 
partition function Z{v,q) for |1 + f| < 1 are all located within a circle given by 
the maximal degree of the graph. The restriction |1 + f| < 1 corresponds to the 
antiferromagnetic Potts model, which is not what we have considered in this paper. 
If the restriction |1 + f | < 1 is relaxed the radius is found to scale as (for spatially 
constant v) 

i?, ~ max [t;, ^ (18) 

where r is the maximum degree of the graph, i.e. the maximum number of edges incident 
on any one vertex. For an ordinary cubic lattice in two dimensions we have r = 4, hence 
we expect to see a crossover from v to v"^ scaling around v = 1. Fig. IHl shows the radius 
Rq as a function of v. 
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Figure 9. The radius Rq for a two dimensional square lattice, i.e. r = 
line is f{v) ~ a • w and the dashed line is g{v) ~ a + 6 • f^. 
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5. Conclusion 

We have shown that the nontrivial information of the Potts model is contained 
in the density rg(C,E), and this is independent of q. rg(C, E) is purely 
combinatorial/geometric property of the underlying lattice, emphasizing the connection 
between these concepts and critical phenomena. Furthermore we have devised an 
algorithm to estimate r£;(C, E) from Monte Carlo simulations, and used this to study 
various properties of the Potts / Random Cluster model. 
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